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Considering Ginsparg-Wilson type fermions dynamically in lattice QCD simulations is a chal¬ 
lenging task. The hope is to be able to approach smaller pion masses and to eventually reach 
physical situations. The price to pay is substantially higher computational costs. Here we discuss 
first results of a dynamical implementation of the so-called Chirally Improved Fermions, a Dirac 
operator that obeys the Ginsparg-Wilson condition approximately. The simulation is for two species 
of mass-degenerate quarks on 12 3 x 24 lattices with spatial size up to 1.55 fm. Implementation of 
the Hybrid Monte-Carlo algorithm and an analysis of the results are presented. 

PACS numbers: 11.15.Ha, 12.38.Gc 


I. INTRODUCTION 

Quantum Chromodynamics (QCD) is a quantum field 
theory with a complicated vacuum structure, involving 
creation and annihilation of particle-antiparticle pairs in 
the vacuum. Presently the most promising approach 
for studying this theory non-perturbatively seems to be 
through Monte-Carlo simulations of the theory on a 
space-time lattice. Fermions, however, are described by 
Grassmann variables and have no simple representation 
appropriate for simulations. Fortunately the fermions oc¬ 
cur bilinearly in the action and in the path integral the 
Grassmann integration can be performed explicitly. This 
leads to the fermion determinant as a weight factor. Each 
fermion species accounts for one such factor. 

In lattice QCD Monte-Carlo simulations the path in¬ 
tegral is approximated by averages over a finite set of 
configurations on a finite euclidian space-time grid. It is 
well understood how to sample gauge field configurations 
with the full dynamics of interacting gauge fields alone. 
It is a challenging task, however, to take into account 
the contribution of the fermion determinant in generat¬ 
ing the probability measure. In the quenched approxima¬ 
tion these determinants are put to unity by hand. Here 
we perform a simulation with full dynamical quarks, i.e., 
taking this determinant into account. 

One of the crucial problems in putting fermions on 
the lattice is preserving their chiral symmetry. As an 
example consider the original Wilson Dirac operator. 
It exhibits spurious zero eigenvalues even at values of 
the fermionic mass parameter that correspond to non¬ 
zero pion masses. This feature complicates the approach 
to the chiral limit. On the other hand, Dirac opera¬ 
tors obeying the Ginsparg-Wilson (GW) condition Q ex¬ 
actly, have no such spurious modes. However, the only 


known explicit and exact realization of such an opera¬ 
tor, the so-called overlap operator Hi is quite expen¬ 
sive to construct and implement in dynamical simulations 
13 Si II Si Various other Dirac operators have 
been suggested, that fulfill the GW condition in an ap¬ 
proximate wav, among them the Domain Wall Fermions 
[ToL ITU IT3 fbl , the Perfect Fermio ns Hi HtH or the Chi¬ 
rally Improved (Cl) Fermions IB Il7| . Such Ginsparg- 
Wilson-type fermions have the chance to allow a better 
approach to the chiral limit, i.e., to come closer to the 
physical value of the pion mass. 

Cl fermions have been extensively tested in quenched 
calculations (see, e.g., Ref. 0 ). In these tests it was 
found that one can go to smaller quark masses with¬ 
out running into the problem of exceptional configura¬ 
tions. On quenched configurations pion masses down to 
280 MeV could be reached on lattices of size 16 3 x 32 (lat¬ 
tice spacing 0.148 fm) and about 340 MeV on 12 3 x 24 
lattices. 

The challenge is to implement these fermions in dy¬ 
namical, full QCD simulations. We discuss here such a 
simulation for two light, mass-degenerate flavors of dy¬ 
namical quarks. For updating the gauge fields we use the 
standard Hybrid Monte-Carlo (HMC) algorithm. We re¬ 
port on our implementation of Dq\ for the HMC updat¬ 
ing and present results of runs at lattice spacings in the 
range 0.11-0.13 fm for 12 3 x 24 lattices of physical spa¬ 
tial extents of 0(1.55 fm). Preliminary results on smaller 
lattices have been presented elsewhere BED]- 


II. ALGORITHMIC CONCERNS 

A. The Dirac operator and the action 

The Cl Dirac operator Dqi was constructed by writing 
a general ansatz for the Dirac operator 


16 

Ai=X4( U)T k , 

k= 1 
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where I\. (k = 1... 16) are the 16 elements of the Clifford 
algebra and a ( U ) are sums of path ordered products of 
links U connecting lattice site i with site j. Inserting D 
into the GW relation and solving the resulting algebraic 
equations yields the Cl Dirac operator Dqi- In principle 
this can be an exact solution, but that would require an 
infinite number of terms. In practice the number of terms 
is finite and the operator is a truncated series solution 
to the Ginsparg-Wilson relation, respecting the lattice 
symmetries, parity, invariance under charge conjugation 
as well as 75 -hermiticity, but connecting sites only over 
a certain distance. For our implementation we use terms 
up to path lengths of four [l|| |T 3 l ■ 

Smearing is an essential quality improving ingredient 
in our Dirac operator. In the quenched studies for Dqi 
UU it was found that smearing the gauge links was 
important (using HYP smearing [2ll l as it resulted in 
better chiral properties for the operator, i.e., the spec¬ 
trum showed less deviation from the Ginsparg-Wilson 
circle compared to the unsmeared case. Therefore we 
decided to use one level of smearing in our studies, too. 
Usual HYP smearing, however, is not well suited for use 
in HMC. The recent introduction of the differentiable 
“stout”-smearing j 2 ^| opened the possibility to imple¬ 
ment HMC for smeared fermionic actions. We thus used 
one level of isotropic stout smearing of the gauge configu¬ 
ration as part of the definition of the Dirac operator. Our 
smearing parameter (« 0.165 in the notation of Ref. ;22|) 
was chosen to maximize the plaquette value. A compar¬ 
ison between HYP and stout smearing on loops of differ¬ 
ent sizes on the same background has been presented in 

El. 


In quenched calculations [la, il7l| it was found that the 
tadpole improved Luscher-Weisz j23| gauge action had 
certain advantages over the Wilson gauge action in the 
sense that the configurations produced with this action 
were smoother than the ones produced by the Wilson 
gauge action. We therefore used this gauge action in our 
dynamical studies. It is given by 




~Pl ^2 q Re tr U P1 ~ q Re tr Ure 

pl re 

-A \ Re tr Utb > ( 2 ) 


where U p \ is the usual plaquette term, U le are Wilson 
loops of rectangular 2 x 1 shape and Utb denote loops 
of length 6 along edges of 3-cubes (“twisted bent” or 
“twisted chair”). The coefficient Pi is the independent 
gauge coupling and the other two coefficients P 2 and P 3 
are determined from tadpole-improved perturbation the¬ 
ory. They have to be calculated self-consistently [IH from 


uo= ( 3 Re tr (kpi) 


1 


a = — 


3.06839 


log (4) 


(3) 


as 


P 2 = (1 + 0.4805 a) , /? 3 = 4 0.03325 a . (4) 

20 Uq Uq 

This determination should be done for each set of param¬ 
eters (Pi ,m). 

Putting together all of the above, the partition function 
that we simulate assumes the form 

Z = J , ( 5 ) 

where M denotes the massive Dirac operator given by 

M(U, m) = Z?ci(m) = D C i(U) + 1 m , ( 6 ) 

U denotes the smeared link variable, and <fi denotes the 
bosonic pseudofermion field. 


B. HMC for Cl fermions 

The alg orithm we used for our simulations is standard 
HMC [2j| as this seems to be the most efficient one for 
simulating fermions at the moment. Although our imple¬ 
mentation follows Ref. [2^| very closely, we believe that 
the additional complications due to the extended struc¬ 
ture of our Dirac operator warrants some discussion. 

HMC consists of a molecular dynamics (MD) evolu¬ 
tion in 2 n dimensional phase space (where n is the orig¬ 
inal dimension of the theory) with simulation time as 
the evolution parameter. This involves introducing con¬ 
jugate momenta and constructing a Hamiltonian for the 
problem. In our case, to construct a Hamiltonian, we in¬ 
troduce traceless hermitian matrices G su( 3) acting 
as momenta conjugate to the (i being the site index 
and /.i the direction of the link) and write 

W = \ tr ( p l2) +Sg + ^ m ) 1 0 ■ (7) 

There are two evolution equations for and 
While Ui'H = , the equation for pi ^ is obtained 

by setting Ti. = 0, 

0 = H = ^2 tr Pj,nPj,v + S g +$ — (M^M) 1 p . (8) 
j, M 

Up to this point there is no difference to HMC with Wil¬ 
son or staggered fermions. The first complication comes 
when we take the time derivative of our operator. Unlike 
Wilson or Staggered Dirac operators, which have only 
one link connecting the neighboring sites, we have longer 
paths. 

As an example let us look at a path of length three. 
Let 




(9) 
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Taking the time derivative, we find 

— i v - JJ. [/t JJ. 

~ 4 All,Ml t7 Jl,Ml U J 2 ,A‘2 t7 J3^3 

+ Uj 2 ,i_i 2 ( — tia)Ujai M3 

(iPh,tiaUj 3l )i 3 ) , ( 10 ) 

where we have replaced f/j,^ by ipi^Ui,^. The time 
derivative of the Dirac operator can thus be expressed 
as 


dM 

dt 


E 


C 3 , M, fc , E ( ^ iPj , M 




( 11 ) 


where IF^ and lW 2 ) contain products of links. To ob¬ 
tain our equation of motion for p.j ^, we have to carefully 
cyclically permute the terms in $ jjj 1 </> until 

all the occur at the same position, at the front of 
the expression. This is of course required due to the 
non-Abelian nature of the variables and it is possible to 
cyclically permute them because we have an overall trace 
over both color and Dirac indices. In fact the trace on 
the Dirac indices can be carried out completely as the 
momenta and gauge part of the Hamiltonian do not have 
Dirac indices. The color trace on the other hand should 
not be carried out as we want an equation for the color 
matrices. 

Finally our equation can be symbolically written in the 
form 


Pj,pPj,n = Pj,ti E ( s ^ a Ples)j,ij, - ipj,n (Dirac). 

( 12 ) 

From this we get our equation for Pj, t i.. To ensure the 
group property for the link variables we actually take 
the traceless part of the r.h.s. of d for Pj tll , as usual. 

Since the D ci contains several hundred path terms, 
it was not practical to do the permutations by hand. 
Thus the major part of the development stage was spent 
in writing automated routines which perform this task 
given the table of coefficients and paths which define the 
Dc i- 

As we have mentioned before, smearing is an impor¬ 
tant part of the Dqi ■ The operator is built not from bare 
links but from smeared links. The momenta, however, 
are conjugate to the bare links. Therefore we have to 
express the smeared finks as functions of bare links and 
compute corresponding derivatives. This is easily done 
for the stout links and we will not discuss it further here 
except to remark that we used a Taylor series with forty 
terms to exponentiate the sum of the staples. While this 
is not the most efficient way of doing the exponentia¬ 
tion, it had the advantage that the differentiation was 
straightforward and the extra overhead in time was neg¬ 
ligible compared to the time required for inverting the 
Dirac operator which is the most time consuming step. 

To sum up some of our operational details, we used 
one set of pseudofermions and the leap-frog integration 


TABLE I: Parameters for the simulations; the first column 
denotes the run, for later reference. The gauge coupling is /3i, 
the bare quark mass parameter am, At the MD time step, 
steps the number of steps for one trajectory. In the last three 
columns we give the acceptance rate in the accept/reject step 
(in equilibrium), the total HMC time of the run and the lattice 
spacing determined via the Sommer parameter. 


# 

L 3 

1 x T 

Pi 

a m 

At 

steps 

acc. 

rate 

HMC 

time 

as [fin] 

a 

~w 

x 24 

5.2 

0.02 

0.008 

120 

0.82(2) 

463 

0.115(6) 

b 

12 3 

x 24 

5.2 

0.03 

0.01 

100 

0.94(2) 

363 

0.125(6) 

c 

12 3 

x 24 

5.3 

0.04 

0.01 

100 

0.93(1) 

438 

0.120(4) 

a 

12 3 

x 24 

5.3 

0.05 

0.01 

100 

0.92(2) 

302 

0.129(1) 

e 

8 3 

x 16 

5.3 

0.05 

0.015 

50 

0.93(1) 

1245 

0.135(3) 

f 

8 3 

x 16 

5.4 

0.05 

0.015 

50 

0.93(2) 

649 

0.114(3) 

g 

8 3 

x 16 

5.4 

0.08 

0.015 

50 

0.94(1) 

776 

0.138(3) 


scheme. In order to speed up the conjugate gradient 
(CG) inverter, we used a chronological inverter by min¬ 
imal residual extrapolation [27| . i.e., an optimal linear 
combination of the twelve previous solutions as the start¬ 
ing solution. This of course was used only in the MD 
evolution and reduced the number of iterations required 
to invert Dqi by a factor between two and three. For 
more details on this see |2Cj| . 

III. SIMULATION 

We tested our code first on 8 3 x 16 lattices m 2(j. 
This lattice size does not allow for small bare quark 
masses, thus we chose am = 0.05 and 0.08. We then 
progressed to the larger lattice size 12 3 x 24 and smaller 
quark masses. Table [I] summarizes the runs discussed 
here. 

The parameters of D ci are fixed by the construction 
principle discussed in Jljj, i.e., by optimizing the set of 
algebraic equations that approximate the Ginsparg Wil¬ 
son equation. In the formalism there is a renormalization 
parameter z s that has to be adjusted such that Dci(m) 
has the low-lying eigenvalues peaked at Re( A) = m. For 
that one studies a sample set of gauge configurations (de¬ 
termined with the full action) and uses that knowledge 
to recursively fix the parameter z s . This recursion can¬ 
not be iterated indefinitely due to limited resources. We 
therefore decided to do this adjustment only at one set of 
parameters (/3i,m) = (5.3,0.05) and than hold all cou¬ 
plings fixed (except for 0 i and m). 

The gauge part of the action also requires a recur¬ 
sive approach as the parameters of the tadpole improved 
Liischer-Weisz gauge action have to be tuned with help 
of the plaquette observable |24|. This can be done by 
equating the Uq to the moving average of the plaquette 
defined over a suitable time interval. The adjustment 
should be stopped once equilibrium is reached so that 
the parameters of the action are fixed. Since we started 
mostly from configurations close to equilibrium, we did 
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not go though this whole procedure, but adjusted the 
value of Uq to be close to its equilibrium value and held 
it fixed. The difference between the assumed plaquette 
and measured plaquette in our runs ranged between 0.2 
and 2%. 

It turned out, that the lattice spacing depends signif¬ 
icantly on both, the gauge coupling and the bare quark 
mass. The bare parameters of the action have no real 
physical significance, however, and we will present most 
of our results in terms of derived, physical quantities. 

The lattice spacing may be determined in various ways, 
e.g., from the values of the measured meson masses or 
from the static potential. At the distance of the Sommer 
parameter ,22], ro = 0.5 frn dynamical quark effects are 
expected to be small and for the masses discussed here 
most likely very small. We therefore use the value of the 
lattice spacing as as derived from the static potential as 
our principal scale (cf. Table [IJ. For the determination 
of the potential we used HYP smeared gauge configura¬ 
tions, since they showed less fluctuation. Details of this 
determination are discussed elsewhere (T§1 |. 

Based on these values for the lattice spacing, our spa¬ 
tial lattice sizes vary from 1.4 fm up to 1.55 frn (for the 
12 3 x 24 lattices) and we may expect finite size effects for 
the derived meson masses. Indeed, meson masses for the 
run (e) (small lattice) are about 10% larger than those 
on the larger lattice of run (d). Also in quenched calcula¬ 
tions at comparable parameters the volume dependence 
of the meson masses was of a similar size when changing 
from 8 3 x 16 to 12 3 x 24 0 . Further increase of the 
lattice size then showed significantly smaller changes in 
the mass values. We also measured the spatial Wilson 
lines but observed no signals for deconfinement. 

The CG tolerance values were fixed to 10~' in the MD 
steps and 10 _1 ° in the accept/reject step. This followed 
from earlier experience on smaller lattices, where we used 
values down to 10 -12 for both (see also the discussion in 
f0). The number of MD steps and the step sizes are 
given in Tabled 

As an example of the equilibration process we show 
in Fig. ^ some observables computed in the run-sequence 
for (3\ = 5.2 and am = 0.03. This run was started from a 
cold gauge configuration and one finds that the plaquette 
quickly approaches its equilibrium value. Another useful, 
more technical, observable is the number of CG iterations 
necessary in the accept/reject step. There equilibrium 
values are obtained somewhat later. 

In Fig. ^ we also plot the development of the lower 
bound to the real parts of the eigenvalues of Dci(0.03). 
The value for every 5th configuration is shown. This 
quantity is a good indicator of the equilibration as it is 
most likely the largest time scale in the problem. At the 
same time it indicates the quality of the Dirac operator. 
Configurations, where the smallest eigenvalue becomes 
zero would lead to spurious zero modes. The amount of 
fluctuation of the observed lower bound also indicates the 
uncertainty width of the choice of the mass. 

In the last row in Fig.^we show the topological charge 
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HMC-trajectories 

FIG. 1: Equilibration signals for (3 = 5.2, am = 0.03; from 
top to bottom: the plaquette mean, the number of CG itera¬ 
tions for the accept/reject step, the lower bound min(Ae(A)) 
for the real part of the eigenvalues of the Dirac operator 
Dfji (rn), the topological charge v. The lower two rows’ values 
have been determined only for every 5th configuration. 

v, again determined for every 5th configuration. The 
definition and determination is discussed in Sect. eh 
S tarting at a cold configuration initially locks the topo¬ 
logical charge to the trivial sector. However, following 
the initial equilibration period, we find satisfactory fluc¬ 
tuation. 

Altogether, based on such observations, we discarded 
the first 100 HMC trajectories before starting the ana¬ 
lyzing measurements. The data are not sufficient to find 
reliable estimates for the autocorrelation length. Judging 
from the inspection of the fluctuation we estimate auto¬ 
correlation times of 10-20. Runs on the smaller lattice 
size had larger autocorrelation times 0]. We analyzed 
every 5th configuration, determining hadron correlators 
and eigenvalues as discussed in Sect. hyi 

All our runs were done with MPI parallelization, typ¬ 
ically on 8 or 16 nodes on the Hitachi SR8000 (at LRZ 
Munich) or on 8 or 16 nodes of two Opteron 248 pro¬ 
cessors at 2.2 GhZ per processor. To get an estimate on 
timing, we note that for the 12 3 x 24 lattice at am = 0.05 
one trajectory takes ~ 4 hours on 8 nodes of the Hitachi; 
one trajectory at am = 0.02 takes ~ 7 hours on 16 of 
the double Opteron nodes. 


IV. WHAT DO WE LEARN 

The numbers that are presented here are affected by 
various limitations: 
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FIG. 2: Plot of the 40 smallest eigenvalues of Dci(0.02) for 
a total of 21 configuration in equilibrium (run (a), lattice size 
12 3 x 24, /Ii =5.2, am = 0.02). 


• We consider only two flavors of mass-degenerate 
quarks. 

• The typical physical lattice in spatial direction is 
only up to 1.55 fm. We thus may expect to see 
finite volume effects. 

• The statistics is limited; most of the runs discussed 
amount to typically 200 to 350 units of HMC time 
corresponding to 40 to 70 independent configura¬ 
tions. 

• The smallest quark mass parameters correspond to 
a pion mass around 500 MeV, the smallest pion 
over rho mass ratio is 0.55. 

• Only a small range of lattice spacing values between 
0.11 fm and 0.13 fm is covered. 

All these constraints are typical for first simulations with 
dynamical quarks, in particular in view of the compu¬ 
tationally demanding properties of the Dirac operator. 
We therefore consider this study as the first but neces¬ 
sary step towards studies on larger and finer lattices with 
better statistics. 


A. Eigenvalues and topology 

Dirac operators obeying the GW condition in its sim¬ 
plest form have a characteristic eigenvalue spectrum: all 
eigenvalues are constrained to lie on a unit circle. The 
operator Dei is only an approximate GW operator and 
its eigenvalues are close to but not exactly on this so- 
called GW circle (radius 1, center at 1 in the complex 
plane). As a first test we studied the low-lying eigenval¬ 
ues of Dei for some gauge configurations obtained with 
dynamical quarks. 



FIG. 3: The distribution of the topological charge for the 
two runs at (3 = 5.2 and m = 0.02 and 0.03. We also show a 
fit to a Gaussian distribution function with the same (a 2 )- 


Fig. |3 exhibits the accumulated eigenvalues for 21 con¬ 
figurations, where we determined the 40 smallest eigen¬ 
values for each configuration. The curve indicates the 
eigenvalue circle for an exact GW operator for mass 
am = 0.02. The bulk of the eigenvalues closely fol¬ 
lows the circular shape. A comparison of the spectra 
on smaller lattices with the quenched case has been pre¬ 
sented in Ref. 0 . 

The complex eigenvalues come in complex conjugate 
pairs Xi and A i and always have vanishing chirality 
(■011 75 1'’/’i)- The real eigenvalues, which would be ex¬ 
act zero modes for exact GW operators, correspond to 
top olog ical charges via the Atiyah-Singer index theorem 
|29l l30l| . For the overlap Dirac operator only zero modes 
either all positive or all negative chirality have been ob¬ 
served. In our case we sometimes (in about 3% of the 
configurations) find zero modes with opposite chiralities 
as well. 

We did compute the chiralities of the real (zero-)modes 
for the runs at /3 = 5.2. The distribution is not symmet¬ 
ric but this is non-significant due to the relatively small 
sample. A symmetrized histogram is in good agreement 
with the Gaussian shape. In Fig. [3 we show the resulting 
histograms together with a Gaussian distribution with 
zero mean and the same second moment ( v 2 ). As is well- 
known from experience with other Dirac operators, when 
evaluated over longer HMC-periods than the ones avail¬ 
able to us, the tunnelmg frequency may show much longer 
correlation time [lj, |3l| ■ 

Tunneling between different topological sectors ap¬ 
pears to be a problem for HMC implementations of the 
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overlap action; various intricate methods have been sug¬ 
gested to deal with it mums As can been seen 
in the lowest row of Fig. U we do not seem to have such 
a problem and observe frequent tunneling. 

As mentioned, we determined the eigenvalues and thus 
the topological charge only for every 5th configuration. 
Comparing the number of changes of the topological sec¬ 
tor we get (for the runs (a,b,d) and (e)) between 12 and 19 
such changes along a HMC-time distance of 100, without 
obvious correlation to the run parameters. E.g., runs (d) 
and (e) have the same values of and m, but different 
lattice sizes. The width of the distribution of the smaller 
lattice is smaller, as expected, but the number of tunnel¬ 
ing events comparable. These values are a lower bound 
for the actual number of tunneling events. In view of 
these we do not find the drastic dependence of the tun¬ 
neling rate on the quark mass observed in simulations 
with the overlap action [?j. 

The topological susceptibility 


Xtop = )/V (13) 

play s a central role in the Witten-Veneziano formula 
13211331 . Due to the axial anomaly the quenched sus¬ 

ceptibility (in the large 7V c -limit) is related to the Ty'-mass 
via 


quenched 

Xtop 



(14) 


(Here f n « 92 MeV denotes the pion decay constant and 
Nf the number of flavors, both defined in full QCD.) 
In quenched simulations values of \top ~ (190 MeV ) 4 
have been found 0, Hi H3 For the dynamical case 
the dependence of ytop on Nf and m has been discussed 
recently, along with a consistent definition of the quantity 
Hi Hi Hi P. Due to the anomalous Ward identity for 
the U(l) axial current the topological susceptibility for 
full QCD should vanish (in the chiral limit). As in formal 
continuum theory j42l 145$ the topological susceptibility 
(now for dynamical fermions) and the chiral condensate 
are related via 


Xtop = + 0(m 2 ) . (15) 

where E denotes the condensate contribution per flavor 
degree of freedom. 

The precision of our results for y t op is too poor to verify 
the linear behavior in the quark mass. For the runs (a) 
and (b) of Tabic [I] we obtain in physical units the values 

Xto P » = (146(8) MeV ) 4 , 

XtoW = ( 166 ( 8 ) MeV ) 4 • ( 16 ) 

In these cases we do find, as expected, that our values 
are definitely smaller than the quenched ones. 


100 


M io 

<D 


1 


0.1 

0 5 10 15 20 

t 

FIG. 4: Correlation function for run (b), for the pseudoscalar 
(top) and the vector (bottom) meson. The curves give the 
results of fits to the cosh-behavior as discussed in the text. 



TABLE II: Meson masses 


# 

aAA 

aM p 

M n /M p 

M n [MeV] 

M p [MeV] 

a 

0.292(10) 

0.535(35) 

0.55(5) 

501(44) 

918(109) 

b 

0.378(8) 

0.619(30) 

0.61(4) 

597(42) 

977(96) 

c 

0.326(18) 

0.502(21) 

0.65(6) 

534(48) 

823(62) 

d 

0.431(8) 

0.626(18) 

0.69(3) 

657(16) 

954(33) 


B. Meson masses 

For the determination of the correlation function we 
used point quark sources and Jacobi-smeared sinks. Our 
definition and notation for the Jacobi smearing followed 
the quenched studies in 0, [45|. We used the narrow 
smearing distribution. Smearing the sink crucially im¬ 
proves the correlator signal quality and is almost as effi¬ 
cient as smearing source and sink. 

The meson interpolating field were the usual ones: 

Pseudoscalar: P = d'ysu , (17) 

A 4 = dj5j4U, (18) 

Vector: 14 = d'fkU, (19) 

where u and d denote the up- and down quark fields; A 4 
is the temporal component of the axial vector, which also 
couples to the pion. 

We computed the correlation functions 

(P(p = 0, t)P(0)) , (20) 

(A 4 (p = 0,f)A 4 (0)) , (21) 

(Vi@=0,t)Vi(p)) . ( 22 ) 

The result (see, e.g., Fig. 0 was then fitted to 

C{t) =D(M) (e~ Mt ±e~ M(T_t) ) . (23) 

The masses thus may be derived from the exponential 
decay and other low energy parameters from its coeffi¬ 
cient. 
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We discuss here only results for the large lattice size 
12 3 x 24. There the fits were done in the range (t a , 4) = 
(6,18). Table El summarizes our results for the masses. 
All error bars have been determined with the jack-knife 
method. We also checked that % 2 /d.o./. for all the fits 
were < 1. The results for the (A4A4) correlator were 
consistent with that for the (PP) correlation function, 
but with slightly larger statistical errors. In the table we 
therefore quote the pion mass obtained from the (PP) 
correlator. 


C. AWI mass and pion decay constant 


The axial Ward identity (AWI) allows one to define the 
renormalized quark mass through the asymptotic behav¬ 
ior of the ratio 


Z A (d t A 4 (p = 0,t)X(0)) 
Zp (P(p = 0,t)X(0)) 


Z m 2m = 2 m^ , 


(24) 


where X is any interpolator coupling to the pion and Za, 
Zp and Z m denote the renormalization factors relating 
the MS-scheme at a scale of 2 GeV. These have been 
calculated for the quenched case at several values of the 
lattice spacing and came out close to 1 [46|. We do not 
know the value for the dynamical case but for the mass 
values presented we expect it to be close to the quenched 
one. We therefore compute the ratio 


(d t A 4 (p = 0,t)P(0)) 
(P(p = 0,t)P( 0)) 


= 2 toawi , 


defining the so-called AWI-mass. 
We measure 


(25) 


(Ai(p = 0,t)P(0)) , (26) 

in order to construct 


(d t A 4 (p = 0 } t)P(0)) . (27) 

Ratios involving the lattice derivative <9 t A 4 depend on 
the way the derivative is taken. Numerical derivatives 
are always based on assumptions on the interpolating 
function. Usual simple 2- or 3-point formulas assume 
polynomials as interpolating functions. We can do bet¬ 
ter by utilizing the information on the expected sinh- 
dependence. In fact, we may use this function for lo¬ 
cal 3-point (t — 1, t, t + 1) interpolation and get the 
derivative at t therefrom. We cannot use correlators like 
(X(t) 9 *^ 4 (0)) since the source is fixed to the time slice 
t = 0 and thus we cannot construct the lattice derivative 
there. 

Eq. (1251) assumes interpolating fields with point quark 
sources. The Jacobi-smearing of the quark sinks intro¬ 
duces a normalization factor relative to point sinks. The 
factors can be obtained from the asymptotic (large t) 
ratios of, e.g., 


_ (P(t)P) _ (M*)P) 

CP (P s (t)P) ’ CA (A^ s (t)P) 


(28) 
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FIG. 5: The ratio 1251 for f3 = 5.2, am = .03. The plateau 
fit range is indicated. 


TABLE III: Results for the AWI-mass in lattice units, and 
the AWI-mass, M n and the pion decay constant in units of 
the Sommer-parameter ro (which is usually assumed to be 0.5 
fm). 


# 

am A wi 

ro Ml 

ro rriAwi 

ro fi r 

a 

0.025(1) 

1.62(28) 

0.103(9) 

0.237(44) 

b 

0.037(1) 

2.29(33) 

0.147(10) 

0.321(45) 

c 

0.037(2) 

1.84(33) 

0.154(13) 

0.314(44) 

d 

0.050(1) 

2.78(14) 

0.195(6) 

0.281(26) 


where the index s denotes the interpolator built from 
smeared sources. We find values of Cp/ca around 0.6 
indicating that the smearing of the quarks affects the 
two operators differently, which may be understood from 
their different Dirac content and “wave function”. We 
did check that the results agree with that derived directly 
from correlators based on point-like quark sinks, only 
the error bars are slightly larger in the latter case. All 
numbers given refer to the normalization for operators 
built from point-like quark sources and sinks. 

Taking this into account we find the AWI-mass from 
plateau values like that shown in Fig.0 The final average 
was taken in the same interval as was used for the mass 
analysis and the error was again computed with the jack¬ 
knife method. The results are given in Table mu 

In full, renormalized QCD the Gell-Mann-Oakes- 
Renner (GMOR) relation relates the pion and quark 
masses: 


/ 2 M 2 = -2mS. (29) 

Here two flavors of mass-degenerate quarks are assumed. 
The quark mass and the condensate (contribution per 
flavor d.o.f.) are renormalization scheme dependent and 
have to be given in, e.g., the MS-scheme. Since the AWI- 
mass is proportional to the renormalized quark mass, this 
linear relationship may hold. Indeed, in lattice calcula¬ 
tions surprisingly linear behavior has been found. 














stant f„. It is related to the coefficient of the (A 4 A 4 ) 
correlator with point source and sink and may be ex¬ 
tracted through its asymptotic (large t) behavior: 

Z 2 a(Aa(p= 0, i)Ai(0)> la ^‘ M n fle- M * 1 . (30) 

Again we assume Za ~ 1; in Fig. 0 we show the re¬ 
sults, which are compatible with the experimental values. 
For this measurement we actually used the point - point 
correlation function. Table 11111 summarizes our results in 
units of the Sommer parameter. 


V. SUMMARY/CONCLUSION 


FIG. 6: The pion mass squared vs the AWI-mass in units 
of the Sommer parameter ro for all 4 data sets (a-d). The 
fitted line corresponds to the lowest order chiral perturbation 
theory behavior, i.e., to the GMOR relation 1 1551 . 
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FIG. 7: f n vs pion mass squared in units of the Sommer 
parameter ro for all 4 data sets (a-d). The * denotes the 
experimental value. 


We have implemented the Chirally Improved Dirac op¬ 
erator, obeying the GW condition to a good approxima¬ 
tion, in a simulation with two species of mass-degenerate 
light quarks. In our simulations on 12 3 x 24-lattices the 
lattice spacing was as small as 0.115 fm. The quark 
masses reached so far correspond to a M n /M p ratio of 
0.55. 

We consider this as a first step towards the final goal 
to approach realistic quark masses. Nice features of the 
Dqi, which were observed in quenched simulations, like 
the spectrum following the GW circle, appear to hold 
in the dynamical case as well. The good consistency of 
the results with experimental numbers seems to indicate, 
that the renormalization factors are close to 1 (as they 
have been found in the quenched case). We also find no 
evidence for spurious zero modes and we therefore do not 
expect serious problems to go down further in the quark 
mass. 

The results are encouraging and demonstrate that the 
Cl operator is well suited for dynamical fermion simula¬ 
tions. 
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